#include <cstdlib>
#include <iostream>
#include <string>
#include <Graph.h>
#include "StabSolvers.h"

// #define NDEBUG
#include <assert.h>

static const char stabstr[] = "STAB";

// initialize variables, constraints and objective function
int Stab::build(StabModeler * bsm, int *erb, int *ere, int *crb, int *cre, int *ub, CPXLPptr solver) {
//	cout << "Starting creating model\n";
//	fflush(NULL);

	Graph * g = bsm->graph();

	// the name of each variable of each subproblem
	int n = g->nverts();
	int m = g->nedges();
	m >>= 1;

//	cout << "stabcols...\n";
//	fflush(NULL);

	int err;
	err = bsm->stabcols();
	assert(!err);
	CPXaddcols(cpxenv, solver, bsm->ccnt(), 0, NULL, NULL, NULL, NULL, bsm->rmatval(), bsm->rhs(), bsm->colname());
	CPXchgctype (cpxenv, solver, bsm->ccnt(), bsm->rmatind(), bsm->type());

//	cout << "edgerows...";
//	fflush(NULL);

	err = bsm->edgerows();
	assert(!err);
	// int CPXaddrows(CPXCENVptr env, CPXLPptr lp, int ccnt, int rcnt, int nzcnt, const double * rhs,
	// const char * sense, const int * rmatbeg, const int * rmatind, const double * rmatval, char ** colname, char ** rowname)
	//
	// ccnt: An integer that specifies the number of new columns in the constraints being added to the constraint matrix.
	// rcnt: An integer that specifies the number of new rows to be added to the constraint matrix.
	// nzcnt: An integer that specifies the number of nonzero constraint coefficients to be added to the constraint matrix.
	// rhs: An array of length rcnt containing the right-hand side term for each constraint to be added to the CPLEX problem object.
	// sense: An array of length rcnt containing the sense of each constraint to be added to the CPLEX problem object.
	// rmatbeg: An array used with rmatind and rmatval to define the rows to be added.
	// rmatind: An array used with rmatbeg and rmatval to define the rows to be added.
	// rmatval: An array used with rmatbeg and rmatind to define the rows to be added.
	//			The nonzero elements of every row must be stored in sequential locations in this array from position rmatbeg[i] to rmatbeg[i+1]-1 (or
	// 			from rmatbeg[i] to nzcnt -1 if i=rcnt-1). Each entry, rmatind[i], specifies the column index of the corresponding coefficient, rmatval[i].
	// colname: An array of length ccnt containing pointers to character strings that represent the names of the new columns added.
	// rowname: An array containing pointers to character strings that represent the names of the new rows, or equivalently, the constraint names.
	*erb = CPXgetnumrows(cpxenv, solver);
	int status = CPXaddrows(cpxenv, solver, 0, bsm->rcnt(), bsm->nzcnt(), bsm->rhs(), bsm->sense(), bsm->rmatbeg(), bsm->rmatind(), bsm->rmatval(), NULL, bsm->rowname());
	assert(!status);

//	cout << "\nisolatedrows...";
//	fflush(NULL);

	err = bsm->isolatedrows();
	assert(!err);
	CPXaddrows(cpxenv, solver, 0, bsm->rcnt(), bsm->nzcnt(), bsm->rhs(), bsm->sense(), bsm->rmatbeg(), bsm->rmatind(), bsm->rmatval(), NULL, bsm->rowname());
	*ere = CPXgetnumrows(cpxenv, solver);

//	cout << "\nedgesymcuts...";
//	fflush(NULL);

	err = bsm->edgesymcuts();
	assert(!err);
//	CPXaddrows(cpxenv, solver, 0, md->rcnt(), md->nzcnt(), md->rhs(), md->sense(), md->rmatbeg(), md->rmatind(), md->rmatval(), NULL, md->rowname());
//	CPXaddlazyconstraints(cpxenv, solver, md->rcnt(), md->nzcnt(), md->rhs(), md->sense(), md->rmatbeg(), md->rmatind(), md->rmatval(), md->rowname());
	if (bsm->rcnt() > 0)
		CPXaddusercuts(cpxenv, solver, bsm->rcnt(), bsm->nzcnt(), bsm->rhs(), bsm->sense(), bsm->rmatbeg(), bsm->rmatind(), bsm->rmatval(), bsm->rowname());

//	cout << "\ncliquerows...";
//	fflush(NULL);

	err = bsm->cliquerows(ub);
	assert(!err);
	*crb = CPXgetnumrows(cpxenv, solver);
	CPXaddrows (cpxenv, solver, 0, bsm->rcnt(), bsm->nzcnt(), bsm->rhs(), bsm->sense(), bsm->rmatbeg(), bsm->rmatind(), bsm->rmatval(), NULL, bsm->rowname());
//	CPXaddusercuts (cpxenv, solver, md->rcnt(), md->nzcnt(), md->rhs(), md->sense(), md->rmatbeg(), md->rmatind(), md->rmatval(), md->rowname());
//	CPXaddsos (cpxenv, solver, md->rcnt(), md->nzcnt(), md->sostype(), md->rmatbeg(), md->rmatind(), md->rmatval(), md->rowname());
	*cre = CPXgetnumrows(cpxenv, solver);

	return m;
}

static double * xperm;
static int xcmp(const void *a, const void * b) {
	double diff = xperm[*((int *)b)] - xperm[*((int *)a)];
	return diff <= 0 ? -1 : 1;
}

static inline void setParam(CPXENVptr cpxenv) {
	// best bound search
	CPXsetdblparam(cpxenv, CPX_PARAM_BTTOL, 0.0);
	// directs messages to screen
	CPXsetintparam(cpxenv, CPX_PARAM_SCRIND, CPX_ON);
	// to avoid warning messages
	CPXsetintparam(cpxenv, CPX_PARAM_MIPSEARCH, CPX_MIPSEARCH_TRADITIONAL);
	// level of printed information
	CPXsetintparam(cpxenv, CPX_PARAM_MIPDISPLAY, 2);
	CPXsetintparam(cpxenv, CPX_PARAM_MIPINTERVAL, 100);
	// set number of threads
	CPXsetintparam(cpxenv, CPX_PARAM_THREADS, 1);
	// used in user cuts
	CPXsetintparam(cpxenv, CPX_PARAM_REDUCE, CPX_PREREDUCE_PRIMALONLY);
	// used in user cuts
	CPXsetintparam(cpxenv, CPX_PARAM_PRELINEAR, CPX_OFF);
	// used in user cuts
	CPXsetintparam(cpxenv, CPX_PARAM_MIPCBREDLP, CPX_OFF);
}

static inline CPXENVptr opencplex() {
	int status;
	CPXENVptr cpxenv = CPXopenCPLEX(&status);
	if (status) {
		printf("CPLEX not open (status = %d)\n",status);
		exit(0);
	}

	return cpxenv;
}

CPXLPptr Stab::opensolver() {
	// base solver
	int status;
	CPXLPptr solver = CPXcreateprob(cpxenv, &status, stabstr);
	printf("CPLEX problem created (status = %d)\n",status);

	return solver;
}

void Stab::setlb(Graph * g, int k) {
	int i, j;
	int n = g->nverts();
	lb = 0;
	if (k == 1) {
		int A[n];
		for (i = n - 1; i >= 0; i--) {
			// assuming symmetric adjacency matrix
			for (j = 0; (j < lb && !g->hasEdge(A[j], i)); j++);
			if (j == lb) {
				S[i] = 1;
				A[lb++] = i;
			}
			else
				S[i] = -1;
		}
	}
	else {
		Graph * gc = g->complement();
		Graphs::dsaturCoverSort(gc, k, S);
		delete gc;
		for (i = 0; i < n; i++)
			if (S[i] != 0) {
				lb++;
				S[i] = 1;
			}
			else
				S[i] = -1;
	}
}

template <>
Stab::Stab(Graph * g, ModelerFactory<StabModeler,Graph> * fm, Weight<double> * w) : StabSolver(), w(w) {
	cpxenv = opencplex();
	solver = opensolver();

	sm = fm->create(g);
	S = new int[g->nverts()];
	setlb(g, 1);

	int ub;
	build(sm, &erowbeg, &erowend, &crowbeg, &crowend, &ub, solver);

	sm->scalerows(erowbeg, erowend-1, 1.0);
	CPXchgrhs (cpxenv, solver, sm->rcnt(), sm->rmatbeg(), sm->rhs());
	sm->scalerows(crowbeg, crowend-1, 1.0);
	CPXchgrhs (cpxenv, solver, sm->rcnt(), sm->rmatbeg(), sm->rhs());
}

template <>
Stab::Stab(Graph * g, ModelerFactory<StabModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix) : StabSolver(), w(w) {
	cpxenv = opencplex();
	solver = opensolver();

	sm = fm->create(g, col, prefix);
	S = new int[g->nverts()];
	setlb(g, 1);

	int ub;
	build(sm, &erowbeg, &erowend, &crowbeg, &crowend, &ub, solver);

	sm->scalerows(erowbeg, erowend-1, 1.0);
	CPXchgrhs (cpxenv, solver, sm->rcnt(), sm->rmatbeg(), sm->rhs());
	sm->scalerows(crowbeg, crowend-1, 1.0);
	CPXchgrhs (cpxenv, solver, sm->rcnt(), sm->rmatbeg(), sm->rhs());
}

Stab::Stab(Graph * g, Weight<double> * w) : StabSolver(), w(w) {
	cpxenv = opencplex();
	solver = opensolver();
	S = new int[g->nverts()];
	setlb(g, 1);
}

Stab::~Stab() {
	delete S;
	CPXcloseCPLEX (&cpxenv);
	if (sm)
		delete sm;
}

int Stab::solveByEnumeration() {
	setParam(cpxenv);
	return CPXmipopt(cpxenv, solver);
}

int Stab::getBest(double * val) {
	return CPXgetobjval (cpxenv, solver, val);
}

int Stab::getNodes(int * nnd) {
	*nnd = CPXgetnodecnt (cpxenv, solver);
	return 0;
}

// MaxStab

template <>
MaxStab::MaxStab(Graph * g, ModelerFactory<StabModeler, Graph> * mf, Weight<double> * w) : Stab(g, mf, w), unit(false), n(g->nverts()) {
	objective();
	sm->delall();
}

template <>
MaxStab::MaxStab(Graph * g, ModelerFactory<StabModeler, Graph> * mf) : Stab(g, mf, new Weight<double>()), unit(true), n(g->nverts()) {
	objective();
	sm->delall();
}

template <>
MaxStab::MaxStab(Graph * g, ModelerFactory<StabModeler, Graph> * mf, Weight<double> * w, int col, const char * prefix) : Stab(g, mf, w, col, prefix), unit(false), n(g->nverts()) {
	objective();
	sm->delall();
}

template <>
MaxStab::MaxStab(Graph * g, ModelerFactory<StabModeler, Graph> * mf, int col, const char * prefix) : Stab(g, mf, new Weight<double>(), col, prefix), unit(true), n(g->nverts()) {
	objective();
	sm->delall();
}

MaxStab::~MaxStab() {
	if (unit)
		delete w;
}

int MaxStab::objective() {
	sm->objective(w);
	CPXchgprobtype (cpxenv, solver, CPXPROB_MILP);
	CPXchgobjsen(cpxenv, solver, CPX_MAX);
	CPXchgobj(cpxenv, solver, n, sm->rmatind(), sm->rmatval());
	CPXchgctype (cpxenv, solver, n, sm->rmatind(), sm->type());

	return 0;
}

int MaxStab::solveByEnumeration() {
	int si = 0;
	sm->startingpoint(S);
	printf("starting point has size %d out of %d expected\n",sm->ccnt(), lb);
	CPXaddmipstarts(cpxenv, solver, 1, sm->ccnt(), &si, sm->rmatind(), sm->rmatval(), NULL, NULL);

	return Stab::solveByEnumeration();
}

// Represent

void RepStab::init(Graph * g) {
	int nreps = rm->nreps();
	repm = new StabModeler *[nreps];
	active = new int[nreps];
	ub = new int[nreps];
	int i;
	for (i = 0; i < nreps; i++) {
		repm[i] = NULL;
		active[i] = i;
	}
	nactive = nreps;
}

RepStab::RepStab(Graph * g, Weight<double> * w) : Stab(g, w), n(g->nverts()) {
}

template <>
RepStab::RepStab(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w) : Stab(g, w), n(g->nverts()) {
	rm = fm->create(g);
	init(g);
}

template <>
RepStab::RepStab(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix) : Stab(g, w), n(g->nverts()) {
	rm = fm->create(g, col, prefix);
	init(g);
}

template <>
RepStab::RepStab(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w) : Stab(g, w), n(g->nverts()) {
	rm = fm->create(g);
	init(g);
}

template <>
RepStab::RepStab(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix) : Stab(g, w), n(g->nverts()) {
	rm = fm->create(g, col, prefix);
	init(g);
}

RepStab::~RepStab() {
	delete ub;
	delete active;
	delete repm;
	delete rm;
}

// Represent

void RepStab::build(int k) {
	int i;

	int ncols = CPXgetnumcols (cpxenv, solver);

	nactive = 0;
	const int * cliquerank = rm->cliquerank();
	const int * rd;
	int nreps = rm->nreps();
	for (i = nreps - 1; i >= 0; i--) {
		repm[i] = rm->submodeler(i, ncols);
		int nr = repm[i]->graph()->nverts();

		if (nr) {
			int j, cr[nr];
			rd = rm->rd(i, &nr);
			for (j = 0; j < nr; j++)
				cr[j] = cliquerank[rd[j]];
			Stab::build(repm[i], &erowbeg, &erowend, &crowbeg, &crowend, &ub[i], solver);
			rm->repcol(i);
			CPXaddcols(cpxenv, solver, 1, 0, NULL, NULL, NULL, NULL, rm->rmatval(), rm->rhs(), rm->colname());
			CPXchgctype (cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->type());
			repm[i]->scalerows(erowbeg, erowend-1, nr, -1.0);
			CPXchgcoeflist (cpxenv, solver, repm[i]->nzcnt(), repm[i]->rmatbeg(), repm[i]->rmatind(), repm[i]->rmatval());
			repm[i]->scalerows(crowbeg, crowend-1, nr, -1.0);
			CPXchgcoeflist (cpxenv, solver, repm[i]->nzcnt(), repm[i]->rmatbeg(), repm[i]->rmatind(), repm[i]->rmatval());
			ncols = CPXgetnumcols (cpxenv, solver);

			if (k == 1 && ub[i] <= lb) {
				int indices = ncols-1;
				char lu = 'U';
				double bd = 0.0;
				CPXtightenbds(cpxenv, solver, 1, &indices, &lu, &bd);
			}
			else
				active[nactive++] = i;
		}
		else {
			rm->repcol(i);
			CPXaddcols(cpxenv, solver, 1, 0, NULL, NULL, NULL, NULL, rm->rmatval(), rm->rhs(), rm->colname());
			CPXchgctype (cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->type());
			ncols = CPXgetnumcols (cpxenv, solver);
			active[nactive++] = i;
		}
	}
	if (k > 1) {
		rm->disjointrows();
		CPXaddrows(cpxenv, solver, 0, rm->rcnt(), rm->nzcnt(), rm->rhs(), rm->sense(), rm->rmatbeg(), rm->rmatind(), rm->rmatval(), NULL, rm->rowname());
	}
	rm->nrepsrow(k);
	CPXaddrows(cpxenv, solver, 0, rm->rcnt(), rm->nzcnt(), rm->rhs(), rm->sense(), rm->rmatbeg(), rm->rmatind(), rm->rmatval(), NULL, rm->rowname());
}

// Max Reps

void MaxRepStab::init(int k) {
	build(k);
	objective();
	int i;
	for (i = 0; i < n; i++)
		repm[i]->delall();
	rm->delall();
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<RepModeler,Graph> * fm, int k) : RepStab(g, fm, new Weight<double>()), unit(true) {
	init(k);
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int k) : RepStab(g, fm, w), unit(false) {
	init(k);
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<RepModeler,Graph> * fm, int col, const char * prefix, int k) : RepStab(g, fm, new Weight<double>(), col, prefix), unit(true) {
	init(k);
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix, int k) : RepStab(g, fm, w, col, prefix), unit(false) {
	init(k);
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, int k) : RepStab(g, fm, new Weight<double>()), unit(true) {
	init(k);
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int k) : RepStab(g, fm, w), unit(false) {
	init(k);
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, int col, const char * prefix, int k) : RepStab(g, fm, new Weight<double>(), col, prefix), unit(true) {
	init(k);
}

template <>
MaxRepStab::MaxRepStab(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix, int k) : RepStab(g, fm, w, col, prefix), unit(false) {
	init(k);
}

MaxRepStab::~MaxRepStab() {
	if (unit)
		delete w;
}

int MaxRepStab::objective() {
	int i;
	for (i = n - 1; i >= 0; i--) {
		if (repm[i]->graph()->nverts() > 0) {
			repm[i]->objective(w);
			CPXchgobj(cpxenv, solver, repm[i]->ccnt(), repm[i]->rmatind(), repm[i]->rmatval());
		}
	}
	rm->objective(w);
	CPXchgprobtype (cpxenv, solver, CPXPROB_MILP);
	CPXchgobjsen(cpxenv, solver, CPX_MAX);
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
	return 0;
}

int MaxRepStab::solveByEnumeration() {
	int i = rm->nreps() - 1;
	int j, nrp;
	while (i >= 0) {
		const int * rp = rm->rp(i, &nrp);
		for (j = 0; j < nrp && S[rp[j]] >= 0; j++);
		if (j < nrp)
			i--;
		else
			break;
	}

	if (i >= 0) {
		int si = 0;
		int nrd;
		const int * rd = rm->rd(i, &nrd);
		int SS[nrd];
		for (j == 0; j < nrd; j++)
			SS[j] = S[rd[j]];

		repm[i]->startingpoint(SS);
		printf("starting point has size %d out of %d expected\n",repm[i]->ccnt(), lb);
		CPXaddmipstarts(cpxenv, solver, 1, repm[i]->ccnt(), &si, repm[i]->rmatind(), repm[i]->rmatval(), NULL, NULL);
	}

	return RepStab::solveByEnumeration();
}

// Rep Colors

class ZWeight : public Weight<double> {
public:
	double weight(const int i) { return 0; };

};

template <>
RepColor::RepColor(Graph * g, ModelerFactory<RepModeler,Graph> * fm, int k) : MaxRepStab(g, fm, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
//	CPXwriteprob (cpxenv, solver, "color.lp", NULL);
}

template <>
RepColor::RepColor(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int k) : MaxRepStab(g, fm, w, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
}

template <>
RepColor::RepColor(Graph * g, ModelerFactory<RepModeler,Graph> * fm, int col, const char * prefix, int k) : MaxRepStab(g, fm, col, prefix, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
}

template <>
RepColor::RepColor(Graph * g, ModelerFactory<RepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix, int k) : MaxRepStab(g, fm, w, col, prefix, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
}

template <>
RepColor::RepColor(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, int k) : MaxRepStab(g, fm, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
//	CPXwriteprob (cpxenv, solver, "color.lp", NULL);
}

template <>
RepColor::RepColor(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int k) : MaxRepStab(g, fm, w, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
}

template <>
RepColor::RepColor(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, int col, const char * prefix, int k) : MaxRepStab(g, fm, col, prefix, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
}

template <>
RepColor::RepColor(Graph * g, ModelerFactory<ARepModeler,Graph> * fm, Weight<double> * w, int col, const char * prefix, int k) : MaxRepStab(g, fm, w, col, prefix, k) {
	rm->objective(new ZWeight());
	CPXchgobj(cpxenv, solver, rm->ccnt(), rm->rmatind(), rm->rmatval());
}

RepColor::~RepColor() {
}

int RepColor::objective() {
	int i;
	for (i = n - 1; i >= 0; i--) {
		if (repm[i]->graph()->nverts() > 0) {
			repm[i]->objective(w);
			CPXchgobj(cpxenv, solver, repm[i]->ccnt(), repm[i]->rmatind(), repm[i]->rmatval());
		}
	}
	CPXchgprobtype (cpxenv, solver, CPXPROB_MILP);
	CPXchgobjsen(cpxenv, solver, CPX_MAX);
	return 0;
}

int RepColor::solveByEnumeration() {
	return 	Stab::solveByEnumeration();
}
